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FOREWORD 


This report consists of a User's Guide and Fortran IV source 
deck for G. D. Sammons' (deceased) Combustion Model Computer Pro- 
gram written under paragraph A.I.3 of Item 0001 in Contract 
F44620-72-C-0046. This report was assembled for publication by 
members of the Technical Staff, Rocketdyne, McGregor, Texas. 


Lloyd R. Lawrence, Jr., Captain USAF, from the Air Force Office 
of Scientific Research, was monitor for the program, 
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SCIENTIFIC REPORT 
MULTIPLE-FLAME COMBUSTION MODEL 
FORTRAN IV COMPUTER PROGRAM 


ABSTRACT 

Several investigators have formulated competing-flame models 
describing composite solid propellant combustion. This report 
covers Mr. G. D. Sanmons' (deceased) competing-flame modeling 
efforts. His papers describing the model are appended. This 
report describes and lists the FORTRAN IV computer program which 
is the practical outcome of any modeling effort. Included is 
the solution of an example problem--estimation of the burning 
rate of a low-smoke HTPB propellant. How the necessary input 
values are obtained and how they are input is discussed in detail. 


Finally, the computer solution printout is depicted. 


The work described here was performed under sponsorship of =? 
the Combustion Energetics Division of the Air Force Office of 
Scientific Research (Contract F44620-72-C-0046). 
Combustion modeling studies will be continued during FY 1975. 
These further studies will be directed primarily toward addition 
of a mathematical treatment for heterogeneous catalysis of burn- 


ing rate to the multiple-flame model. 
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INTRODUCTION 


Considerable efforts over the past seventeen years have been 
directed toward elucidating the mechanism of combustion of com- 
posite solid propeilants. The ultimate goal of these efforts has 
been a truly comprehensive and descriptive mathematical combustion 
model, The most sophisticated models stemming from these efforts 


(1) (2) 


(3)(4) 
Derr and Price (BDP) . These combustion models allow one to 


are the competing-flame models of Hermance and Beckstead, 

use some independently determined kinetie and thermal parameters 

to predict burning rate. This report covers one of the more recent 
modeling efforts--Mr. G. D. Sammons' (deceased) extension of the 


BDP competing-flame model. 


Mr. Sammons! papers presented at the 9th and 10th JANNAF Com- 
bustion Meetings, which describe his model, are appended to this 
report as Appendix I. His computer program originally written in 
FORTRAN IV for use at a time-share computer terminal has been re- 
written for use on an IBM-300. This rewritten program is described 


here briefly. Also given here is a solved sample problem, 


COMPUTER PROGRAM 


The NRCOMB Computer Program, listed completely in Appendix B, 
is written in FORTRAN IV. [t has been checked out on an IBM 500 


and has been used to estimate the burning rate of composite solid 


propellants having carboxyl-terminated and hydroxyl-terminated poly- 


butadiene binders. 


A flow-chart of NRCOMB is reproduced in Figure 1. 
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FIGURE 1 


SUBROUTINE FLOW CHART OF j 
COMBUSTION MODEL COMPUTER PROGRAM 
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NRCOMB is the main program. Its sole function is to control program 


flow. Calculations are carried out in the subroutines, 


Subroutine INPUT does as its name implies--inputs the parameters 
used to calculate burning rate. These parameters are listed in 
detail in the section "SAMPLE PROBLEM: PROGRAM INPUTS USED TO 
COMPUTE BURNING RATE." 


Subroutines QUTPUT, ENERGY and FACTOR are output routines that 
print out parameter names and values and column heads burning 


rate vs. pressure data, flame parameter vs. pressure data. ele, 


Subroutine CONCAL increments the pressure from a start value 
to a stop value, 

Function subroutines BOB and ROB establish the functional 
relationships XNAME DON 1 and XNAME = SURFAC, 

Function subroutines CONVR 2 and CONVR 3 are convergence rou- 
tines that force convergence on the desired value when upper and 
lower extremes are known, CONVR 2 yields surface temperature; 
CONVR 3 yields critical AP particle size. 

The group of subroutines and function subroutines CONVR, ROB, 
CONVR 3, CRIL, OXID and PART utilize the multi-modal oxidizer 
particle size input to compute an average particle size, critical 
AP size, fraction of AP greater and less than critical, and an 
oxidizer ignition time. The roles of ROB and CONVR 3 have been 
described. CRTL computes a critical AP size; @XID computes an 
average diameter and a greater and less than critical size weight 
fraction, PART computes the AP ignition time utilizing the AP 


particle size distribution. 


Subroutine STAMP 1, utilizing several subroutines for the 


competing flame, computes surface temperature and as a final 


RO Oe or he oe ae 


—_ 
Re eee Te 


i 
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step calculates burning rate using the expression: Total Mass 


Flow/Propellant Density. STAMP 1 calls on BOB and CONVR 2, both 


: 
i 
4 
| 


of whose roles have been described. It also calls on DON 1 which 


in turn calls on SOXCAL, NEWRAF, FLAME, BESSEL and BROOT. 


Function subroutine DON 1 computes propellant surface tempera- 
ture. It uses S@XCAL to compute the fraction of total prepellant 
surface involved in AP only decomposition, It utilizes FLAME, 
BESSEL, BROOT and NEWRAF to compute the parameters needed in 


primary diffusion flame and final flame calculations. 


SAMPLE PROBLEM: PROGRAM INPUTS USED 
TO COMPUTE BURNING RATE 


Given here are descriptions of the program inputs used to 
estimate propellant burning rate. For this example, a typical 
low-smoke HTPB propellant having the composition listed below has 


been selected. 


Ingredient Weight Percent 
HTPB Binder 12.0 
Ammonium Perchlorate 
400 / (nominal diameter) 55.5 
200 & (nominal diameter) 13.0 
80 - (nominal diameter) 12.2 
7-11 & (finely ground) 5.0 
Aluminum, 5 __ 0.3 
100.0 
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Oxidizer Particle Size and Weight Fractions - 
INPUT 


S_DMAX, DMIN, AND SIG 


The oxidizer particle size distribution given in the formula- 
tion is broken down into several particle size weight fraction 
ranges for inputing by using tne typical AP particle size analyti- 
cal data given in Table 1. Such an input for the example pro- 


pellant is shown here. 


DMAX (KL) DMIN(/) Subscript 
600 4209 (8) 


420 9306) 
14,06) 


gol) 
604) 
yg) 


5(e) 


From nominal 400% fraction 

From nominal 400P and 200 fractions 

From nominal 200 and 802 fractions 

From nominal 200K, 80 and 7-ll- fractions 


From nominal 80% and 7-11) fractions 


The 600-420 fraction is obtained by multiplying 0.555 by 0.37; 
the 420-230 fraction by multiplying 0.555 by 0.59 and 0.13 by 0.23 
and using the sum, etc.; to produce a 7-modal AP distribution 

(NN = 7) input. 
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TABLE I 


TYPICAL AP PARTICLE SIZE DISTRIBUTIONS 


A. Nominal 400 4 AP (as received) 


Particle Diameter Range 


Microns Weight Percent 
600-420 37.0 
420-350 29.0 
350-250 30.0 
250-210 1.0 
210-177 1.0 
177-74 2.0 


B. Nominal 200 4 AP (as received) 


Particle Diameter Range 


(Microns) Weight Percent 
500-230 23.0 
230-140 59.0 
140-78 13.0 
78-25 5.0 


C. Nominal 80 4 AP (as received) 


Particle Diameter Range 


Microns Weight Percent 
110-72 60 
72-54 20 
54-35 10 
35-11 8 
11-4.3 2 


D. Fine Grind - Rocketdyne's Nominal 7-11 4 AP 


Particle Diameter Range 


Microns Weight Percent 
80-60 5 
60-9.8 75 
9.8-5 15 
Dag y) 
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Modality of AP - Input NN 


In the previous section, the oxidizer was divided into a 


7-modal distribution. So this input is: NN = 7 


Weight Fraction Oxidizer and 
Aluminum - Inputs ALFA and EPS 


From the formulation, it is evident that the 
values are: ALFA = .857 
EPS = .003 


Density of Oxidizer, Binder and 
Aluminum - Inputs RHOX, RHOF and RHAL 
Values for these quantities are (standard values): 
RHOX = 1.95 gm/ew? 
RHOF = 0.91 on /ene 
RHAL = 2.70 gm/cm? 


Initial Propellant Temperature and Flame 
Temperature - Inputs TZERO and TF 


Flame temperature value is obtained from thermodynamic calcu- 
lations. Values fo. these inputs for the example propellant are: 


TZER# = 300 K 
TF 2908 K 


Binder and Oxidizer Kinetic Parameters - 
Inputs AF, EF, A@OX and EPOX 
Pyrolysis rates for the separate fuel and oxidizer components 


are expressed by Arrhenius expressions of the form: 


th = Aexp (E/RT.) 
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where A and E are input parameters. For the fuel component, the 
9 

frequency factor AF = 299 em/cm™ sec and activation energy 

EF = 16,900 cal/mole which were obtained by pyrolyzing an HTPB 


(7) 


binder are used as inputs. For the oxidizer component, 


ammonium perchlorate, A@X values ranging from 104 to oe" and 
EGX values ranging from 17,000 to 73,400 for pure AP have been 
reported. For use here, we have selected frequently used values 
in modeling studies--namely, AX = 4 x 10° gm/cm> sec and 


EX = 30,000 cal/mole. 


Gaseous Diffusivity - D, Input as GAMMA 


In this program, gaseous diffusivity is calculated using 


J 
Beckstead, et al's axeseien © i) 


1.75 
T 
y= D, E P 


Do is an input parameter. It is estimated using the Chapman-Enskog 
(5) 


relation : 


Dip = 0.0018583 


and 


P 

° z a <a 

Po PAB ¢ a | os Pe 
0 
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For this sample calculation, we have assumed that the gaseous 


fuel species is C which is consistent with the PMW value given 


318 
in a later section. For oxidizer we have used 0,. Using appro- 
priate values for these gases in the Chapman-Enskog equation, a 
value for C,H,-0, diffusion of 0.0928 ou /aee is obtained which 


38 2 
yields the value of dD, below. 


D, = 5:07 x 107° (input as GAMMA) ch’ /aee 


Dap might also be estimated using the empirical relationship 
developed by Gilliland which can be found in many chemical engi- 


neering textbooks 


Ammonium Perchlorate Oxidizer 
Flame Temperature — Input TAP 

This is the temperature of one of the three competing flames-- 
the temperature of the oxidizer monopropellant flame. The standard 


(3) (4) 


value used by Beckstead is used here. 


TAP = 1400 K 


Stoichiometric Ratios of the Primary Flame 
and of the Final Flame - Inputs XNUP and XNU1 


To estimate these stoichiometric ratios, the general formula 


below is used: 


(Stoichiometric Ratio) = odes ue) Suet ae 
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AP decomposition can be represented by the stoichiometric 


equation: 


NH, C10), —> ‘0.5 N,, + 0.8 HCl + 1.6 H,,0 Oe Cl Ae De 


‘) 


-~ 


which gives a molar oxygen flow rate of Tly.5 


= 0.0102 moles 
oxyzen/gm 
The approximate formula for an HMDI cured Poly BD R-45M HTPB 
binder (the formula used for material balance in thermodynamic 
calculations) provides a basis for estimating the moles of fuel 
per gram of fuel flow. Using this formula and setting the heaviest 


species at CL gives the following: 
( 


3 Co Species 
(mole w/~ 300)* <> 7 C, Species 


H N 
CoL.i5 31.95°.753 258 
OS 16 CH,, Species 


--about 8.7 moles of fuel species per 300 grams of binder on the average so 


that 
ANUP = Dita oy 
0.01 
* Polybutadiene is {CH, 5x 
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In the final diffusion flame further breakdown of the binder decomposi- 


tion products is assumed so that 
aap 7 Cs Species 


91 15451.95°.753" 258 9%, 16 CH, Species 


This is, on the average. 11.35 moles of fuel species. This gives 
fa J 


— 11.5/300 
UL ee = 9.85 


The diffusion flame solution employed in this model requires XNUP and XNU1 


values greater than 2.8. 


Molucular Weights for the Primary and Final Flames - 
Inputs PMW and GMW 


To be consistent with stoichiometric ratios these values are calculated 


00 + 1l 
rv " = 
Mole Wt. Moles Fuel + Moles 0, 


as follows: 


giving 


} ue 
GMW Oe; © 33 (Average mole wt. of combustion 
ie products of this propellant 
is approx. 30) 


Burning Surface Temperature--Input TS 


The program calculates the temperature of the burning surface. This 


eclculation requires an initial estimate. For use in this example: 


TS = 750 Kk 
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Flame Height Factor and Oxidizer Ignition Dela 
Time Parameters--Inputs AFH, CIGN, POWIGN, and PWD 


AFH 0.5 

’ mM, on 
CIGN = 190 sec atm /cm 
PQWIGN = 0.75 (an exponent: m) 


i 


POWD 0.8 (an exponent: n) 


The flame height factor and ignition delay time parameters are 


(1)(2)(3)(4). 


those used by both Hermance and Beckstead 


Rate Constant Parameters for the Diffusion 


and AP Flames--Inputs KPF, KAP1l and KAP2 


These are the rate constants for gas phase reactions of 
arbitrary reaction order and expressed in the following general 


form 


Values used here for these k's are 


KPF = 30 gm/cm? sec Atm 
KAP] = 1.12 gim/cm sec Atm 
KAP2 = 2.13 gm/cm? sec Atm 


These are "standard" values for the gas phase rate constants. 


Reaction Order Constants for Diffusion 


Flame and AP Flame--Inputs XNI1, XN2, and XN3 


XN1 = 1.5 
XN2 = 1.8 
XN3Z = 1.8 


S 


These values should be considered "standard" values for the 


parameters in AP composite propellants. 
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Average Ce Value for Flames and Solids and 
perceived | SDE RAE ee ra ee ee aeD 
Gas Conducti.ity - Inputs CSUBP and XLAMB 

me is used in this model in several calculations. For 
example: 

(1) The energy balance term accounting for propellant heat up 
mC (Eat) 
(2) As an enthalpy term in the AP flame relationship 
C(t p-T,) 

(3) In non-dimensional standoff distance calculations in 

the multipliers 

C CM 

e He and _p ox 
p- oe 
kach requires a a value. This model, like the BDP model uses 


only one C, value for both gases and solids. This value is an 


overall average value. Why this approximate approach has been 


used becomes evident on calculating this overall average value. 


The specific heat of AP beiow and above the phase transition 


temperature yields an average value 
( 
c= 800 + 585 | 9.537 cal. /em-K 


The specific heat of the example propellant is 


cr = 0.352 cal/gm-K 
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Each gaseous product contributes ics proportional share to or 
To obtain the value for the AP flame, the products shown in the 


stoichiometry section provided the basis for the materia) balance 
and the AP flame temperature and propellant surface temperature pro- 


vided the average temperature from which the eS contributions in the 


table were obtained. 


2 
T(ave) = Lo ee = 1075 k 


2 


Gas Wt. Fraction cy Contribution” 
Ny 0.096 0.0259 

HC] 0.201 0.0386 

H,0 0.198 0.1026 

Cl, 0.241 0.0300 

05 0.264 0.066 


* Estimated from charts in Hougen and Watson, "Chemical 
Process Principles: Part One - Material and Energy 
Balances," John Wiley 


These contributions sum up to the ae value for the AP flame. 


C = 0.263 cal/gm - K 
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To cbtain the oy value for the propellant flame, the final 


combustion product distribution was used to previde an estimated 


average value at an average temperature as shown below: 


faves \< Surface Temp - Flame Temp. & 750 : 2910 | 1830 K 


Gas Wt. Fraction* a Contribution** 
No 0.102 0.0287 
H,0 0.286 0.1604 
C0, O47 3 0.0500 
co 0.174 0.0491 
HC1 0.265 0.0560 


* From flame temperature calculations 


** Estimated from charts in Hougen and Watson, "Chemical 
Process Principles: Part One - Material and Energy 
Balances," John Wiley. 


These tabulated contributions sum up to 


c. = 0.344 cal/gm - K 


15/R-4827 


AN Rockwell International 


In arriving at these “A values for the flames, ideal gas 
behavior has been assumed. 


The overall average of these values is used as input. 


2 + 0,263 + 0.3! 3 ; 
csupp = 2-352 + 0.55 . Sey Os. 0 ge) cal/gm > k 


Beckstead's "Standard" mi? for gas conductivity (ideal 


behavior assumed) is used here as gas conductivity input. 


XLAMB = 0.0003 cal/cm sec - K 


Aluminum Heat Absorption - Input QAL 


The latent heat of fusion - the heat absorbed when molten 


drops of aluminum are formed on the propellant surface--is input. 


QAL = [(Heat of Fusion) + €, (mean) x oo] = 257 cal/gm 
Heat_to Vaporize the Fuel 
Binder - Input QFUEL 
(7) 


For this input the measured value for an HTPB binder is 


used. 


QFUEL = 433 cal/gm 


Heat Stemming from AP Decomposition, Reaction of Less Than Critical 
Size AP with Fuel, and from Liquid Layer--Inputs QL, QSOLID and QLIQ 
The BDP model employed a QL input only to account for decompo- 


(3) (4) 


QL, was varied to account for more or less reaction in the solid 


sition of the oxidizer, AP, at the surface This parameter, 


phase. Mr. Sammons introduced two additional adjustable parameters, 
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QSOLID and QLIQ, to account for reaction in the solid phase. Here 
the following values for these parameters (values used by Sammons) 
are input. 

QL = 0.0 cal/gm 


QSPLID = 1085.0 cal/gm 
QLIQ 1000.0 cal/gm 


u 


Pressure Programming Inputs PSTART and PSTOP 


Two pressure programming inputs to start and stop the incre- 
mental pressure increases employed in calculating a burning rate 


vs. pressure curve are required. Inputs used here are: 


PSTART = 1.0 
PSTOP = 220.0 


lt 


Pressures used in the burning rate calculations are in atmos~ 


pheres. These pressures are converted to psia before outputing. 


SAMPLE PROGRAM RUN 


Following are the computer inputs and computer printout 


obtained with the example data described in the previous section. 
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CARD 1: TITIE CARD 

NO. REQUIRED: One card per run 
FUNCTION : Identify run output 
FORMAT : (20A4 ) 


79 columns of alphanumeric data may 
be input to identify run 
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CARD 2 
NO. REQUI 
FUNCTION: 
FORMAT: 


sw bt ¥ tert 


nternationa 
CONTROL CARD 
RED: One card per run 


Specify run options to be executed 


(1216) 


Specifies run option 
1 - Input complete set of data 
(NJPB must always be 1 for first run) 
2 - Stacked run; input new values for 
AF and AQx 
- Stacked run; 
- Stacked run; 


3 input new value for QL 
4 input new value for XNU1] 
5 - Stacked run; input new value for AFH 
6 - Stacked run; input new values for 
PSTART and PSTUP 
7 - Stacked run; input new values for 
IX and TZEKO 
8 - Stacked run; input new value for QF 
9 - Stacked run; input new modality of 
oxidizer and aluminum weight fraction 
10 - Rerun same problem - not used 
11 - Rerun using weight fraction and 
modality of oxidizer data stored in 
INPUT1 
100 - Repeat printout of last run 
105 - Repeat last run - printout only HDN1, 
0--Do not write tape for plotting 


HDP1 and S@X] 
7--12 IPLAT If IP1p 
input 
2--Write log(P) and log (rate) on logical 


13--18 16 IPLT 
unit 10 for plotting. Do not print 
40 


output. 
fy ) “Write log (P) and log (rate ) on 


T > 0, do not print propellant data 


logical unit 10 for plotting and 
print output 


Not read - may be used to identify card 
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CARD 3: DATA CARD 

NO, REQUIRED: One card per run for NJOB=1 
FUNCTION: Specify propellant parameters 
FORMAT : (6£12.5) 


TZERO Initial MRE temperature , 
deg K 


| aura | Oxidizer weight fraction 


ie ee 


Flame temperature, deg K 


Molecular weight of final flame 


Stoichiometric ratio of the final 
flame > 2.8 


Mean AP particle size, 
0.0 for multimodal AP 


Not read - may be used to identify 
card 


yw input as 
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CARD 3A: DATA CARD 


NO. REQUIRED: One card per run for NJOB=2 through 8 
(Do not include for NJOB=1 or NJOB > 8) 


FUNCTION : Specify input data changes to previous run 
FORMAT : As noted 


E12.5 ieee cae Frequency factor--See Card 4 NJOBe2 
13--24 E12.5 ca Frequency factor--See Card 5 


]--12 E12.5 sey from AP decomposition, NJOB=3 
cal/ gm 

]--12 E12.5 XNU1 Stoichiometric ratio of final | NJOB=4 
flame > 2.8 


E12.5 | AFH | Flame height factor NJOB=5 


]--12 E12.5 PSTART Pressure to start incremen- 
tal calculations, atm 
NJOB=6 


BE12.5 PSTOP Pressure to stop incremen- 
tal calculations, atm 
ee er ie: as Output control(l 
se Se 


13--24 E12.5 TZERO Initial propellant tempera- 
ture, deg K 

1--12 E12.5 QFUEL Heat to vaporize the fuel NJOB=8 
binder, cal/gm 


(1) 1 oN 


ormal cutput 
2 Print energy terms vs pressure in sub- 
routine energy 
IX = 43. Print weight factors with respect to 
critical particle size and critical par- 
ticle size vs pressure in subroutine factor 
Print burning rate and XTAR's vs pressure only 


= 
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GY Bocestsyre Dovtein 
Rockwell international 


CARD 4: DATA CARD 

NO. REQUIRED: One card per run for NJOB=] 
FUNCTION: Specify propellant parameters 
FORMAT : (6E12.5) 


]--12 E12.5 QFUEL may) to vaporize the fuel binder, 
cal/ ar aa 


13--24 E12.5 RHOF Density of binder, | Density of binder, go/em) 


Frequency factor, (i team -sec 


57-48 | Activation energy, (!)cal/mole 


49--60 Stoichiometric ratio of the final 
flame > 2.8 
61--72 E12.5 | pw | Molecular weight of primary flame Molecular weight of primary flame 


73--80 Not read - may be used to identify 
card 


(60 ayaa eer 


inetic parameters used to calculate binder pyrolysis 
rates were: 
= AF exp(EF/RTs) 
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OR ecctore nc 


WW CARD 5: DATA CARD 
NO. REQUIRED: One card per run for NJOB=1 
FUNCTION: Specify propellant parameters 
FORMAT : (6E12.5) 


Density of oxidizer, gm/ cm 


Frequency factor(1), gm/cm?-sec 


Activation bierey lt): cal/mole 


Ammonium perchlorate flame tempera- 
ture, deg K 


Not read - may be used to identify 
card 


Miinetic parameters used to calculate oxidizer pyrolysis 
rates were: 


m = AOX exp(E0X/RTs ) 
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\ 


. 
13--24 
37--48 E12.5 PSTART Pressure to start incremental cal- 

culations, atm 


Rocketdyne Division 
Hoar eebattiac 


CARD 6: DATA CARD 

NO. REQUIRED: One care per run for NJOB=1 

FUNCTION: Specify oxidizer ignition parameters and pressure 
range for calculations 

FORMAT : (6F12.5) 


Oxidizer Gantereienn cacy hes time co- 
efficient\), sec-atmPOWIGN/ 
cmPOWD 


POWIGN Oxidizer er tty delay time pres- 
sure exponent 


Oxidizer Lava ays delay time diam- 
eter exponent 


PSTOP Pressure to stop incremental cal- 
culations, atm 


ae 


Not read - may be used to identify 
card 


ee 


xidizer ignition delay time parameters use in the following 
expression: 


beds fg fp tek | Gee” 
2 pPOWIGN 
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CARD 7 


ON Fecctre nsnn 


DATA CARD 


NO. REQUIRED: One card per run for NJOB=1 


FUNCTION : 


FORMAT : 


Specify rate and reaction order constants 
for diffusion flame and AP flame 


(6E£12.5) 


Colum 


Rate coefficient for diffusion 
flame, gm/cm3~sec-atm 


Rate coefficient #1 for AP flame, 
gm/cm>-sec-atm 


Rate coefficient #2 for AP flame, 
gm/cm3-sec-a tm 


Reaction order constant for diffu- 
sion flame 


Reaction order constant #1 for AP 
flame 

Reaction order constant #2 for AP 
flame 

Not read - may he used to identify 
card 
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ON Fesessre ono 
Rockwell International 


CARD 8: DATA CARD 

NO. REQUIRED: One card per run for NJOB=1 

FUNCTION: Specify thermal properties and flame height 
FORMAT : (6£12.5) 


Pees [om Gas conductivity, cal/cm-sec- -deg K 


Gaseous diffusivity, em2/sec 


Flame height factor 


73--80 Not read - may be used to identify 
card 
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ca a 


IR fesetsre nse 


Lf CARD 9: DATA CARD 
NO. REQUIRED: One card per run NJOB=1 
FUNCTION: Specify propellant parameters 
FORMAT: (6E12.5) 
rales eg korea 
]--12 E12. Initial estimate of cee uae surface 
temperature, deg K 
13--24 E12.5 QSPLID Heat from reaction of less than 
critical size AP with fuel, cal/gm 
25--36 E12.5 QLIQ Heat from liquid layer, cal/gm 
37--48 E12.5 Aluminum heat absorption, cal/gm 
49--60 . Density of aluminun, gm/cm? 
73--80 Not read - may be used to identify 
card 
{ 
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OI cctore noes 


CARD 10: DATA CARD 

NO. REQUIRED: One card per run for NJOB=1 or NJOB=9 

FUNCTION: Specify modality of AP particle size distribution 

FORMAT : (1216) 


Description 


Modality of AP particle size dis- 
tribution 


Not read - may be used to identify 
card 
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/ 


eee 


CARD 11: DATA CARD 


NO. REQUIRED: One card for each two particle distribution 
modes input for NJOB=] or NJOB=9 


FUNCTION: Describe AP particle size distribution 
FORMAT: (6E12.5) 


1--12 E12.5 SIG(1) Weight fraction of AP particles be- 

tween DMIN(1) and DMAX(1) 

13--24 E12.5 DMIN(1) Minimum diameter of particles in 
this mode, microns 

25--36 E12.5 DMAX(1 ) 

37--48 E12.5 S1G(2) Weight fraction of AP particles be- 
tween DMIN(2) and DMAX(2) 

49--60 E12.5 DMIN(2) Minimum diameter of particles in 
this mode, microns 

oe oe San 


Maximum diameter of particles in 
this mode, microns 


Maximum diameter of particles in 
this mode, microns 


Not read - may be used to identify 
card 
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md 


stn nic eae 


OR rece ncor 


CARD 12." DATA CARD 

NO. REQUIRED: One card per run for NJOB=1 or NJOB=9 
FUNCTION: Specify weight fraction aluminum 
FORMAT : (6£12.5) 


13--72 E12.5 eee, Not used 


73--80 Not read - may he used to identify 
card 
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Qe = 


CARD 13: DATA CARD 

NO. REQUIRED: One card per run except for NJOB=105 
FUNCTION: Specify print option 

FORMAT : A4 


1--4 A4 IAN Specifies printout option--YES 
blank in columns 1--4 implies 
printout only rate vs pressure. 
Any other entry will cause print- 
out to include surface temperature, 
CRIT, n/Do, M@X, SPX, and XTAR vs 
pressure to be printed. 
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OR mr | 
| 


CARD 14: TITLE CARD 

NO. REQUIRED: One card per run 

FUNCTION: Specify end of job or title of next run 
FORMAT : (20A4 ) 


79 columns (2--80) of alphanumeric 
data may be input to identify 


stacked run. 


END OF JOB in columns 1--10 implies 
end of job and computer will ter- 
minate normally. 


NOTE: If stacked job is to be run, repeat input indicated for Cards 2 
through 13 as required by NJOB value specified in Card 2. 
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RESULTS OBTAINED WITH EXAMPLE DATA 


- 7 MONNAL AP = 86% SOLIDS = 90.37% ALUMINIIM 


BURNING RATES 


PROPELLANT DATA IS 


273.8 MICRON 


- w 
° ic: 
1! oa coal { 
c o Ww 
° o 

io} - 

i Oo 
" Vs 

( 2833 e 


iC hm ee ell @ 
oumer a. 
ee Zs 


b> LL Lt CD ee 


OQ 
e 
Cc 
[7 
— o 
wil 
< i aa 
- ee 
<q oad 
Cc Ww 
wz 
oO 
z ze 
-_ G2 
4 wee 
[a 4 ze 
=} CG 
c aw 
Cc 
?¢ 
Cf 
z 
G 
ert 
a 
| nl 
4 
Ge 
las 
a 
uw 
~! 
_ oa 
Cc ee 
— oo 
=< fog 
GS ~ 
tou 
ZO 
ox 
om CT 
wo 
33/R-4827 


ACTIVATION FNFRGIEFS AND RATE FACTNRS ARF 


Cu 
eouun 
oo 


cNWO oO 
—~ eO UT 
Ones 
aca 
Wee 
aco 
UuUaara 
waCcCc 


FLAME PROPEPTIES ARE 


0.5900 


0.00039 


1400.0 


27.0 DEG CENTIGRADF 


a) 
hee 
a 
z 
Li 
- 
ad 
< 
a 
— 
pot 
Zz 
~—_ 
— 
z 
< 
a 
- 
ul 
a 
a 
a 
a 


Colt 


SUPF T 


a 


rc 


Vi 


IN/ 


RQ, 


CAS SFC 


ac, 


OSTA 


PT, 


PT ATMS 


OKC AHA ROR NK MHCO 
TNE OR SEC AAR MIN! 
SNR OACOG CIOL CH 


EYIOSTAR MET VTAOUR OST 
Cat lT oT POMS COT ORD 
eeeeeee#es8ee#eseset 8 
HAR LESTORO—MGIN A OCLC 
Cm ete RAG OOeh we 
INU UNE IN IFIn oor Oo 
NO Cer or arceororre 
eocrwncacnu avon 
LUE CHMNC OU HSI OAC 
COCO tet eT COST 
@eee#ee##v#2e#e%s8etee 8 @ @ 
Sceenaneaccoao0o0cecen 
ACEC erucoroo ce uce 
CROCCAFROCIU NRO Sea 
PE CTOMMOC KOR ANS Se 
RM CINIAIE PUM HORA 


x) 
e e 
OCCOOCOCOOCO COO Heirs 


cCcoaccoccoccccooc 
eeoeeee7#nree#ese ® 8 e888 8 & @ 
MOcOcoocococoaooeoe 
Onw CwWOCoCCECCCOOCO 
SKM AOMILEM CAL CUCO 

Saat IAI AIM SY 


CC CMHC OS UE O ms 
AGHINVOOTA SHOLCOM a 
ee@e#eeertesesteeeeese 
Va CMmroOor sRS RMN COLN 
MK HVUAIM SLA CME OPM 
mM, 


HNHe HDN “rx sn 


OT ATM 


WWE OM OOFAON 
MAM PINE COE EE 
RRARAR RR RR AR 
MOraceravr eae MmoOrmn 
@eeoe0o0e#e@e##8#@¢e¢e 
oO 


OOoeracoo0o0oe 


Occacoooeoco0o0oco 
eccoocooococooo0oea 
ocococc¢coooecoo0es 
Val alValtallal tal aital'al'alallaltal allel’ s 
@eesvoeve@eeeeve7ee 8 Gee @ 
ocooocoococooeeecoo 
prrrprerdreroreenne 


osccececcaceeocoo 
ialetelealelelaialelalelolelelole) 
eoccococcoccoococo 
Val ea Stel ai al Vel aly al’ al lal! allel lalla’ a) 
eececceoesteseoeoeeese 
ocoocce 
tt t 4 


OCD etree A OSIM © a 
Far OOLROGOCLOACC] 


*eeeoeeee2eoeoesuesecesve#eree8 @ee 
GACCKCE IF CaNOOCIN 
SMa ANAM SC OOME OR 


Katain ci 
34/R-4827 


XTAR PF XTAR AP XTAR D 


XTAR PD 


B®, IN/SEC 


PSIA 


PTs 


SOK DL BMMCOKMMIE OC O 
LAO CALE AIEEE ON em aN 
coeoeveeee eee ee eee 
RANNMMS SUEOFZOOON 
ot et 


BANE OneCare einecg OE NW 
CAMO COLTOAEACKHME EO 
eeeeeeensteee e088 80 
SCHPROHCNONS (NO 
PINS MAKER 


RWC IMWOCEOR SOHCLS 
LAS COP LTOR CUS OH Ne 
eeooeereeroeo eee ea eee 


WNLONNN HH HOCTCOOOC 


TORE NINOSRE SOO 
NDete SMM MN MOO OMNT 
eeeoeevoteeveeveeeeveeoee @ 
OQOLOONK Mh LKATMOP 
STHEKKWOAHMSEGMONOGE 


NOOHCrereceraqocseor 
COORUNCT CONWAY Oa 
COR C HMOC Mets a san 
COCO Stet enim cre 
@eeeeses8nhkreeee#s8e@ 


a 
maococaccccocoocs 


ft 
« 
‘Oo 
“ 


0.5 


qooooco0oeooocococoe 
eeeeodoeeeeseeeveeesr 
oooooacococooca0o0 
onwmer acoocoococe 
HAHANNALOPR OAWOMOCS 

SMtNN + 


35/R-4BP7 


Divielon 


7 ere ine 


Combustion modeling studies will be continued during FY 1975. 


These further studies will be directed primarily toward addition of 


a mathematical treatment for heterogeneous catalysis of burning rate 


to the multiple-flame model. 


(1) 


(2) 


(4) 


(5) 
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ADVANCEMENT OF COMPREHENS IVE 
COMBUSTION MODELING 


G. D. Sammons 
Member of Technical Staff 
Rocketdyne Solid Rocket Division 
McGregor, Texas 76657 


INTRODUCTION 


Studies in many laboratories over the past 15 years on the mechanism 
of combustion of composite solid propellants have provided enough informa- 
tion that it has been possible recently to approach the goal of a truly 
comprehensive and practical mathematical combustion model. This paper is 
a progress report covering the first five months of effort under AFOSR 
Contract F44620-72-C-0046 with this goal as its aim. 


Atte the pioneering effort of Prof. C. E. Hermance,! Beckstead, Derr 
and Price“ advanced this type of modeling by introduction of their multiple 
flame model. The purpose of the current work is to advance another step 
forward in comprehensive modeling. 


In both of the earlier models, researchers found that a relatively 
large amount of condensed phase or surface reaction was necessary. This 
is entirely consistent with a large volume of experimental evidence gath- 
ered in recent years.?’*'?»”% A concept of critical particle size for ammo- 
nium perchlorate (AP) was presented recently? that allows the actual cal- 
culation of condensed phase heat release. This theory is based on the 
concept that a critical particle size exists for AP such that all parti- 
cles beneath this size are consumed in condensed phase reaction. 


MATHEMATICAL TECHNIQUE 


The critical particle size is arrived at by calculating at a specific 
pressure the thermal wave penetration into the solid phase sufficient to 
cause a very fast condensed phase reaction. Initially in this study, a 
simple limiting thermal gradient was calculated by dividing the burn rate 
by the propellant thermal diffusivity and multiplying by the difference 
between surface temperature and a temperature at which the rate of con- 
densed phase reaction becomes rapid enough to occur in the residence time 
available. 


Figure 1 shows the surface energy balance incorporating the critical 
particle size model. The symbols Moxy and Sox are defined as the mass flux 
and burning surface of supercritical oxidizer. Mp is the mass flux of fuel 
binder plus subcritical oxidizer. Delta is the weight fraction of this 
interstitial propellant that is consumed in condensed phase reaction and 
is calculated based on stoichiometry considerations. 


The flame heights are calculated by the Burke-Schuman treatment as 
modified for short flames by Beckstead.- A planned improvement in this 
calculation is to use more terms of the series solution. Use of the first 
term gives only a rough approximation. / 


Equations actually used in the computer program are presented in 
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Pig. 2. The first i8 syeacy to Raise ENERGY TO VAPORIZE ENERGY TO VAPORIZE ENERGY FROM 
an equation for sur- PROPELLANT TO T, OXIDIZER FUEL CONDENSED PHASE 
face temperature that 

is derived from the Son s .) 
surface energy bal- ™%y°>%! * - ™ 3 % - 20) ™ z Sie bn, z a 


ance using mass bal- g 


ance equations. The 
coefficient gamma is ENERGY FROM PRIMARY ENERGY FROM FINAL 
the weight fraction DIFFUSION FLAME ENERGY FROM AP FLAME DIFFUSION FLAME 
of supercritical 4 5 5 

oes ETM ‘s B 
crassa in the prem «fog a (Sy) 0-800 [nally owe) 
puter program iter- 


ates on surface Figure 1. Critical Particle Size Model Energy Balance 
temperature with an 

estimated value of crit- SURFACE TEMPERATURE 

ical particle size and ieee en Se au - 6a os iy er Oi i, : 
after convergence a Foo Boek k Cy t Cs 

OE eae ~~ is cal- B, Op, 

culated. e surface exp (-€,) + 1 - Bo) ex ( a? 
temperature is then re- ; Sor ‘ % i ae bud 
calculated if the ini- oan Ce’) 

tial and calculated oS L 

critical particle sizes PRIMARY FLAME 

are different. This Qe OT Te YQ + A-OK - 6a -VIQ, 


doubly iterative proc- 


ess is continued until FINAL FLAME 


apie C 
convergence of critical Soe [x Soh eevee 1] u-6)-Y¥), . 6u-¥) 
particle size is ob- FY ts hoe aciap alee Y e tes 
tained within 0.01 
micrometers. The equa- Figure 2. Equations Used in Computer Program 


tions for the unit heat 

release in the primary and final flame are also shown in Fig. 2. They also 
include the calculation of condensed phase heat release and the coefficient 
delta, 


COMPUTER PROGRAM 


There are several reasons that a time share computer was used for this 
work. One is speed of turnover; for instance, 15 or 20 turnovers are ob- 
tained within 30 minutes as opposed to one or two a day with conventional 
hardware. It allows direct communication with a large, fast computer and 
the editing software is extremely versatile. Additionally, it provides a 
way for an engineer not versed in computer technology to easily run a con- 
plex program with varying inputs. It is hoped that this program will 
eventually be a tool for engineers and chemists. 


Figure 3 shows two types of input that are most commonly used. There 
are many other options] inputs what can be used if desired. The input on 
the left uses a data file called FLEXI to input about 30 ingredient paran- 
eters. For this option a description of oxidizer must be typed in, which 
includes modality and weight fraction and minimum and maximum particle size 
for each mode. At this point either output can be asked for or any single 
parameter can be changed; for instance, in this example the unit heat 
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release in the condensed phase was changed. Only the information after 
each question mark is typed,as input all the rest is output by the computer. 


RUN 
NRCOMB 14853 NR T/S JULY 21, 1972 
NJOBs IPLOT, IPLT 


? 12000 


TITLE LENGTH. IDENTIFY THE RUN 
? O2TEST 


1S DATA FILE TO BE USED 2? YES 
FILE NAME ? FLEX! 
MODALITY OF OXIDIZER 7? § 


SIGCJ) 2» DMINCJ) » DMAX( J) NKCOMA OR:47 NR T/S AUGe O09. 1972 
? «212230, 450, 


? «3521502210, 

? 0122750142, NJOBsIPLOT, IPLT 
2? ©2562 1504,6202, ? 1161.0 

2? ©0603+8s68+9 


FILE NAME ? FLEXI 


NJWOBs IPLOTs IPLT 
? 4,020 NJOB, IPLOTs IPLT 
? 1001.0 


QSOLID= ? 750 


NJOB, I[PLOT.s IPLT 
7? 1022.0 


Figure 3. Computer Printout of Theoretical Data for CTPB Propellant 


The input on the right of Fig. 3 is a minimum input where all data are 
entered by the computer. Again, iz any one parameter change is desired it 
can be input as was done with QSOLID in 


the other input. This is the input RATES ONLY 7? YES 
that is used for the study of param- PTs PSIA BR _IN/SEC CRIT 
etric variations. 
100-0 001799 10-000 
120-0 0-192 ° 

The simplest form of output is 150-0 acaGsl bt 
shown in Fig. 4. Only burn rate pres- be ect ie 
sure and critical particle size are 300.0 0-8369 a 
given. Particle size is in micro- 400-0 0.2521 8.168 
meters. If no subcritical AP is cal- che prea ee 
culated, the critical particle size 1000-0 023336 60465 
column does not appear. Additional ces @ pai 5-955 
data can be obtained without any recal- engeen oleae par 
culations simply by inputting the dad polit ab 
proper index. Some of the tables of $000.0 ockeas Saabs 
parameters are shown in Fig. 5 through NJOBs IPLOTs [PLT 
8. Figure 5 lists the coefficients Trees, 
used in the surface temperature and 1X= 73 
other equations shown in Fig. 2. The ; 
parameters SIGC and SIGSC are the ba gers el eal 
weight fractions of subcritical and 
supercritical AP, respectively. The Figure 4. Continuation of Printout 
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sum of these is always one. 
tions to surface temperature in degrees Kelvin. 
vaporization are negative contributions. 


in micrometers. 


the total time for about twice the output given in Fig. 5 through 8. 


DELTA 


02266 
02261 
0.252 
02243 
02236 
0236 
0-228 
0-222 
00207 
00176 
OelS% 
06127 
0-081 
00041 
0-005 
0-000 


Figure 6 and 7 show the calculated contribu- 
The oxidizer and fuel 


Figure 8 gives the flame heights 
The computer time of 30.5 seconds shown in this figure is 


SIGC 


0-060 
0-058 
10056 
0-053 
0-052 
0.082 
0-049 
0-048 
0-044 
0-036 
0-031 
0-025 
0-015 
0-007 
0-001 
0-000 


S1GSC 


0+«940 
0.942 
0.944 
0°947 
0-948 
0-948 
0.951 
0-952 
0-956 
0-964 
0969 
e975 
0-985 
0.993 
0-999 
1-000 


CAMMA 


0-790 
0.791 
02793 
06795 
0-797 
0.797 
0+ 799 
0-800 
0-803 
0-810 
Oe614 
O-B19 
06827 
0-834 
0-839 
0-840 


Continuation of Printout 


IN/SEC JXTD VAP 


Di» PSEA BHe IN/SEC  JeTAF 
100+0 OelI9) 160000 
12060 061929 140000 
15060 022041 160000 
200-0 002247 160000 
250-0 0+2353 1+0000 
300-0 062369  0+8090 
400-0 062521 05028 
5000 022603 063444 
700-0 22822 062007 

1000-0 0+3330 061209 
1200-0 0+3710 00952 
1500.0 0+4269 0+0734 
2000-0 065246 0+0548 
2500-0 Oe6171 060456 
3000-0 0°7066 0+0407 
4000-0 0-8685 00247 
VOR, LPLOTs [PLT 
? 7990 
1x3 22 
NJOBs IPLOT, IPL 
? 100.020 
Figure 5. 
PT» PSIA Bue 
100-9 001799 
120-0 0.1929 
15060 0+2081 
200-0 0.2247 
250-0 0-2353 
306-0 0+2369 
400-0 0-2521 
500-0 0+2603 
700-0 0+2832 
1000-0 0+3336 
1200-0 0-3710 
1500-0 0-4289 
2000-0 0-5246 
2500-0 0+6171 
3000-0 0+7066 
4000-0 0+8685 
Figure 6. 


Figure 9 shows the sigmoid distribution curve usually obtained with 


commercial oxidizer. From this curve it is easy to see that the use of 50 


0-0 
0-0 
0-0 
0-0 
0-0 
0-0 
0-0 
0-0 
0-0 
0-0 
0-0 
0-0 
0-0 
000 
000 
OeN 


FUEL VAP 


0.2 
Oc2 
O«2 
0-2 
O52 
0-2 
Ne2 
O«2 
Oe2 
O52 
Oe2 
0-2 
0.2 
0-2 
02 
002 


CP HEAT 


10409 
102.2 
978 
9302 
9001 
9001 
8509 
8303 
7601 
6269 
54-8 
42-9 
2602 
12.7 
eS 
0-0 


Continuation of Printout 


CRIT 


19-000 
90635 
90208 
&.789 
8-522 
8.522 
Be1l68 
10948 
70395 
60465 
5-955 
5289 
42457 
30887 
32470 
2-899 


CRIT 


10-000 
92635 
9208 
8-789 
6.522 
8.522 
8-166 
710948 
72395 
6.465 
5-955 
5-280 
4.457 
3.687 
3-470 
2-899 


percent point does not give a true representation of particle size. 
improve realism in oxidizer input, a spline fit was made to the sigmoid 
curve and the data input as a multimodal oxidizer. 


was input as four modes. 


single mode if necessary. 
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To 


The curve in Fig. 9 


The oxidizer handling subroutine separates the 
modes into supercritical and subcritical particle sizes and will divide a 


PIs PSIA Blts IN“ SEC PF HEAT APF HEAT FF HEAT 
100-0 0.1799 43209 0-0 Oe0 
120-0 001929 44104 seas) 0.0 
150-0 00-2081 452-5 0-0 0-0 
200-0 00e2247 464-5 0-0 0-0 
250-0 02353 47205 0-0 Oe) 
300-0 0+2369 47202 203 Ool 
400-0 0.2521 461-0 2202 008 
500-0 0-2603 41904 6709 1.8 
700-0 0+2832 31609 162.0 206 

1000-0 063336 21864 305¢1 1.0 
1200-0 03710 17903 359-4 0-3 
1500-0 0-42289 143-4 416064 Ool 
2000-0 0+5246 Lilel 478-5 0-0 
2500-0 06171 9404 519-6 0-0 
3000-0 0+7066 85-5 549-0 000 
4000-0 00-8685 5301 597-9 0-0 
NJOBs IPLOT. I PLT 
? 72000 
IX ms | 
NJOBs IPLOTs IPLT 
? 1002020 
Figure 7. Continuation of Printout 
PTs PSIA BRs IN/SEC NTAR PD XTAR PF XTAR AP 
100-00 0-1799 11-79 14616 205-99 
120-00 0-1929 12-41 11-54 162-44 
150-00 0-208! 13-01 8-91 120-59 
200-00 062247 13-76 6-25 80657 
250-00 02353 14¢25 468 58-00 
300-00 002249 14059 3059 42093 
400-00 0.2521 15-02 2048 27266 
500«01 02603 Sell 1663 19618 
700-00 0-283? 15206 1-20 11.28 
000-00 003336 14-86 0-83 6-82 
200-00 00-3710 14-72 Oe 70 5°37 
500-00 0+4289 14-26 0-58 4-07 
2000-01 025246 13-25 0-46 2-88 
2500-03 066171 12-08 0-39 2-22 
3000-01 0+ 7066 10-83 0.34 1-81 
4000-02 08685 12-40 0-27 1629 
NJOBs IPLOTs IPLT 
tS 
RUNNING TIMEs 30-5 SECS 1/40 TIME 3 18-4 SECS 
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BYE 


OFF AT 13250 


Figure 8. 


Continuation of Printout 


RESULTS 


CRIT 


10-000 
9-635 
9-208 
8.789 
8.522 
6-522 
8.168 
72948 
70395 
6-465 
5-955 
5-260 
4-457 
30887 
3470 
2.899 


XTAR D 


2553 
27222 
29-39 
32012 
34-06 
34-82 
37016 
38-29 
40-6 
46-27 
§0-25 
56.47 
66-87 
717016 
87234 
104.33 


The initial test of the computer program was made with the data of 
Bastress primarily because both Hermance and Beckstead used these data. 
It, therefore, provides a point of correlation between the models (Ref Fig. 
10). The theoretical data from the critical particle size model are shown 


as points. Data for the 20-micron AP (monomodal, 70% loading) fits Bastrees' 


data within experimental error. The fit is not as good with the 200-micron 


AP but_it is a considerable improvement over Beckstead's last published re- 
sults.2 One of the major factors making the 200-micron monomodal propellant 
a better fit was to use the entire oxidizer distribution curve as input, as 
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was described previ- r— 
ously. No further work  !®} 
is planned with this 90 | 
propellant, as it is 
not one currently in 
wide use. All subse- 
quent comparisons of 
theoretical data will 
be with production or 
development propellants, 


see 


GROUND AP 15 


Figure 11 shows 
the theoretical data 
obtained for an 84% 
AP-loaded propellant 
Wath or CaEboyy % 2 3 456 8 10 2 3 456 810 2 3 456 
terminated poly- DIAMETER IN. MICRONS 
butadiene (CTPB) binder. 

Experimental data show Figure 9. Particle Size Distribution of Ground 
are for an 86% AP load- and of Commercial 200-micron AP 

ing. This line would 

ae Oe eps ” 8 EXPERIMENTAL (BASTRESS] 


at high pressure if —-~ THEORETICAL (BECKSTEAD) « 7u 
there were only an Oo} THEORETICAL 20 1G |e Gh a 


84% loading. The fit THEORETICAL 200 
of data is, there- 05 T ; 
fore, much better 
than it appears to be 
in Fig. ll. 


% WEIGHT LESS THAN DIAMETER 


_ ———— 


tN. SEC 
o 
= 


RATE, 


All experimental 
data in Fig. 10 and 
1l are strand rates, 
which are not neces- 
sarily predictive of 
motor rates. In the aan | 
future only motor vt 
rates will be used 0 ee 
for correlation if 200 mm 500 700 1000 2000 3000 4000 
at all possible. PRESSBRE,: Ets 


BURNING 
J 
mn 
4. 


Figure 10. Correlation of Theoretical Burn Rate 
with Bastress Data for Polysulfide 
Propellants 


CONCLUS IONS 


Work completed thus far indicates that a comprehensive and predictive 
steady-state combustion model of composite solid propellants is possible. 
These models are complex but modern computing techniques make these models 
easy to use and study. It is our philosophy in this study to be realistic 
and include as much known phenomenology of solid propellant combustion as 
possible hut to avoid unnecessary complexity. It is also our purpose to 
program this model so that it can become a tool for propellant chemists 
and engineers. 
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The next step in 0.8 , SEES SE Ba GO GO — 
this program is to add ~~ — EXPERIMENTAL r+ 4 ht “+ 
aluminum, The equa- Oot 7 THEORETICAL i; —_ Ss 4 4 


tions have already 


been written and the et i f 1 ¥ T 7 
treatment will be much o 5 4, 2 | | 
the same as is used ny ? 
by Cohen, Derr and z 
Price. es ata if T 
= 

Currently, gas 2 
phase conductivities & 0? + Se 
are calculated using = | | 
Chapman-Enscog 
theories and then 
input to the program. | 
Later these calcula- | | am Piet) ey Pa | , L 1 
tions will be written m0 300 500 700 1000 2000 3000 4000 


into the program. PRESSURE, PSIA 


Figure 11. Theoretical Burning Rates for an 84 
The flame height Percent Oxidizer Loaded Propellant 
calculations will be improved with CTPB Binder 
by using more terms of the 
series solution. 


Parametric studies run continuously with this program are indicating 


that the exothermic liquid layer on AP during burning is very important. 


A treatment should be added based on the surface area of the burning super- 


critical oxidizer. This surface area is calculated witnin the program. 


The effect of catalyst is to be added by changes in activation energies 


for the various chemical reactions treated by this model. 
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NOMENCLATURE 
- surface temperature 


- initial propellant temperature 


- weight fraction of supercritical AP particles in total 
propellant 


- specific heat of gas phase 
- latent heat of vaporization of AP 


~- weight fraction of binder/critical AP particles consumed 
in condensed phase reaction 


- energy of fuel decomposition 
- energy of condensed phase heat release 
- fraction of reactants reacting in primary flame 


- nondimnensional flame standoff distances for monopropel lant 
flame, primary flame and total oxidizer flame 
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ABSTRACT 


This progress report covers the first 15 months of a continuing study 
on construction of a comprehensive mathematical model of composite solid 
propellant combustion, with emphasis on work during the past 12 months 
Since reporting on the early efforts at the meeting last year. 


The liquid layer believed to exist on burning ammonium perchlorate 
has been treated in a manner similar to that described last year. Depth 
of penetration of a thermal wave is calculated and this value is multiplied 
by the surface area of oxidizer on the burning propellant surface. The 
surface is calculated in the computer program for use in other equations. 
Total volume of the liquid layer derived in this manner is then used to 
calculate the thermal contribution to surface temperature. 


Treatment of aluminum has been added, and the treatment of multimodal 
oxidizer has been improved. Most internal parameters have been thoroughly 
checked to verify that they have realistic values. Considerable time was 
spent on the Burke-Schumann flame height calculations. It has been found 
that at least seven terms of the series solution must be considered. The 
Single first-term calculations that have been used in the past do not give 
a good assessment of flame height; the flame height is expressed implicitly 
in the final equation, and it was necessary to solve it by a Newton-Raphson 
iteration. Improvement in a number of input parameters is described. Gas 
phase conductivities were calculated using Chapman-Fnskog theories with 
the Hirschfelder-Eucken corrector. The same theories were used to calculate 
diffusion coefficients. Stoichiometric ratios were calculated using more 
realistic pyrolysis species. A thermogravimetric analysis was made of a 
polybutadiene propellant binder and a statistical method was used to derive 
kinetic parameters from the thermograms. The method used was the Zsake 
modification of Doyle's method, which is probably the most sophisticated 
way of treating TGA data. Some other improvements in programming and cor- 
relation with actual propellant motor firings are shown and future plans 
are given. 


INTRODUCTION 


Work continued during the past year on the comprehensive modeling of 
composite solid propellants.* The model being developed is based on the 
Beckstead, Derr, Price multiflame model (1,2) and utilizes a critical par- 
ticle size theory (3) to provide a means of calculating condensed phase 
heat release. A number of improvements have been made to the model and the 
computer program that was described at the meeting last year (4); one of 


“this work is being supported by AFOSR under Contract F44620-72-C-0046. 
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the major improvements being the addition of a treatment of the liquid 
layer on burning ammonium perchlorate (AP). These improvements, an exan- 
ple experiment and future plans are described in the following paragraphs. 


MATHEMATICAL MODEL AND IMPROVEMENTS 


The most recent surface energy balance is shown on Fig. 1. The new 
terms (not shown last year) are the energy to melt the aluminum and the 
energy from the liquid layer. In this equation m,, is the mass flux of 
interstitial propellant; or, more specifically, the binder plus the sub- 
critical AP. The term for condensed phase heat release based on the. 
critical particle size theory is included. The critical particle size 
theory, stated simply, is that there exists for a given propellant and 
pressure a particle size where all smaller particles will be consumed in 
condensed phase reaction. The terms for aluminum fusion and liquid layer 
on AP are also present. The treatment of sluwinum used is somewhat dif- 
ferent from that used in the BDP model; with the methodology explained 
later. 


ENERGY TO RAISE ENERGY TO VAPORIZE ENERGY TO VAPORIZE ENERGY FROM ENERGY 
PROPELLANT TO Yr OXIDIZER FUEL CONDENSED PHASE TO MELT AL 


$ s 
. 1-6im £ Q + 6m, = 


mc (T, * 1) . -m -—=— Q 4 


ENERGY FROM PRIMARY ENERGY FROM FINAL ENERGY FROM AP 
DIFFUSION FLAME ENERGY FROM AP FLAME DIFFUSION FLAME LIQUID LAYER 


* Be Opp me (é, Seg a3 * be wo(E,) a ah é, ) * My ‘ 19 


Fig. 1. Critical Particle Size Model Energy Balance 


Figure 2 shows first the current equation for surface temperature. 
This was derived from the energy equation using mass balance equations in 
the same manner as originally used by Beckstead (2). The following three 
equations are for the heat release in the primary, final, and AP flames, 
as denoted by the subscripts (PF, FF, AP). 


Flame Heights. The Burke-Schumann equation modified for short flames is 
shown in Fig. 3. In the past a single term of this expression has been 
used but we have found that at least seven terms should be used. The equa- 
tion is really not as complicated as it first appears; it is nothing more 
than is shown in the second equation. By adding the extra terms the equa- 
tion has been made implicit in flame height but it is solved easily by a 
Newton-Raphson iteration. Convergence is obtained rapidly when the initial 
guess for flame height is taken as one-tenth the value obtained from a 
single-term calculation. This table is the numerical value of the first 
five terms; and each term is about half the value of the preceding term. 

An actual check of flame heights when different numbers of terms are used 
shows that at least seven terms wust be used to get the flame within 1 
micron. 
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Fig. 2. Current Equations for Surface Temperature 
Used in Computer Program 


FLAME HEIGHTS 


l 
v-(sy) ce ye ze J (Cad (i eva,e]’ <a} 3 
[, 


exp j- 
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A * By exp (-C)N) + B, exp (-CoM) +... By exp (CN) 


SIZE_OF TERMS 


1 - 0.072 
2 - 0.086 
3 - 0.024 
4- 0.011 
5 - 0.005 


SOLVED BY NEWTON-RAPHSON ITERATION 


Fig. 3. Burke-Schumann Equation Modified 
for Short Flames 
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A number of other improvements have been made to the computer pro- 
gram for this model. During the past year all parameters have been 
thoroughly checked, particularly those used internally and normally not 
listed. We found that some of the volume and weight fractions would be- 
come unrealistic; e.g., greater than one. When this occurred it caused 
discontinuities in the burn rate curve. This was remedied by setting 
boundaries on these parameters. 


Treatment of the oxidizer particle size has been improved. The spline 
fit of the particle size distribution curve is still used. A normal bi- or 
tri-modal oxidizer blend will usually be split into 10 modes by this tech- 
nique. Calculation of the h/D ratios is handled in a different manner 
than that described by Cohen, Derr, and Price (5) last year (Fig. 4). Each 
mode of supercritical oxidizer is split into 50 parts; the ignition delay 
term in the equation for h/D is calculated for each part, summed and 
averaged. Each ignition term is then multiplied by the weight fraction 
of the mode and summed with similar terms produced from each mode. Then 
bh/D and Sox are calculated. 


SURFACE DIMENSIONALITY 


n 

r KD 
Beet sabe) fi - za He 
) 7A \(1 " * Fox pm. 


DIFFUSION COEFFICIENT 


LIQUID LAYER TREATMENT 


s. 0. 2 
. oo we Mx 
oe Te, 


Fig. 4. Calculations for h/D Ratios 


Oxidizer/fuel stoichiometric ratios are different in the primary flame 
and the final flame. At temperatures such as are found at the binder sur- 
face it is known that polybutadiene depolymerizes (6) with the major specie 
produced being butadiene. This is considered the major diffusing specie 
in the primary flame, giving a stoichiometric ratio of 3.7. The Chapman- 
Enskog equation in the center of Fig. 4 was then used to calculate the dif- 
fusion coefficient for the primary flame. In the final flame it was 
assumed that hydrocarbon fragments were the major diffusing species; a 
stoichiometric ratio of 8.0 was used, along with a larger diffusion co- 
efficient. \ 


Liquid Layer. The AP liquid layer was not treated as a melt but as an 
exothermic layer. We have always assumed this layer to be a solution of 
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AP and perchloric acid. A thermal wave penetration is calculated in the 
Same manner as for the critical particle size. This volume is converted 

to weight fraction of AP in the propellant that is available for exothermic 
surface reaction (Ref equation at bottom of Fig. 4). 


EXPERIMENTAL RESULTS 


The combustion model with the improvements noted above was applied to 
several smokeless propellants for which experimental data were available. 
Theoretical calculations for three of the propellants (based on carboxy- 
terminated polybutadiene binders) were made at two initial propellant 


temperatures and temperature sensitivities of burning rates were calculated. 


As shown in Table I, agreement of 77 F burn rates is good. Slopes of the 
burn rate curves are identical. The SO, values were calculated using the 
one over one minus "n" relationship. The theoretical 9, values are all low 


and should be two to three times higher; and we will look into this further. 


Teble I. Temperature Sensitivity 


EXPERIMENTAL 
r at 1000 psi, 


THEORETICAL 
ket at 1000 psi 


in./sec aes %/deg F in./sec Sere %/deg F 
0.8 0.38 #60.10 0.16 0.38 0.36 0.04 0.06 
0.29 0.50 ‘0. 7 0.33 0.27 0.50 0.02 0.05 
0, 45 05 O10 0.21 0.4 0.50 0.03 «80.07 
Figure 5 shows a 0.8 
burning rate vs pres- 
sure plot; the upper 0.6 


line represents a least 

squares fit to rocket 

motor firings; the lower eo 
line, the theoretical. a 
This is a two cycle -by- z 
two cycle log-log plot. ws 
The scales are the same = 
as those widely used in z 
industry to plot rocket 3 
motor rate; therefore, 

the data are presented 

in a realistic, prac- 

tical manner. As noted 

last year we did not use 0.1 
strand data for corre- 
lative purposes since 
strand rates are not 
necesSarily predictive 
of motor rates. Pro- 
pellant used in the motor from which the plotted data were taken has a 
Butarez CTL-type binder and an 86% solids loading. The trimoda] oxidizer 
is comprised of commercial] 200- and 80-micron fractions and Ma jac-ground 


PRESSURE, PSIA 


Fig. 5. Propellant with Butarez CTL Type with 


86% Solids and Trimodal Oxidizer 
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10-micron AP at weight percent levels 35/26/25. The theoretical curve is 
identica] in slope up to 1000 psi. 


Comparable data 0.8 
with a propellant con- 
taining a trimodal 
mixture of 400-, 200-, 
and 80-micron AP at 0.5 
weight percent levels 
51/20/15, respectively, 
are shown in Fig. 6. 
The break at 1000 psi 0.3 
is well verified by 
motor firings. No 
curve fitting was used 
to obtain these re- 
sults; i.e., they were 
obtained on the first 
computer run. Fig. 7 
shows data from a 1 
hydroxy-terminated 200 300 500 700 1000 2000 3000 4000 
polybutadiene (HTPB)- PRESSURE, PSIA 
based propellant with : 

a large amount of fine Fig. 6. Propellant with Butarez CTL Type 
oxidizer; i.e., a tri- Binder with 86% Solids and 

modal blend of 200-, Trimodal Oxidizer 

80-, and 10-micron at 
weight percent levels 
of 31/15/40, respec- 
tively. Again, we are 
comparing motor data; 
and, although the burn 
rate 18 somewhat lower, 
the two lines have 
identical slopes. Data 
from a propellant that 
has a tetramodal oxi- 
dizer, i.e., 400-, 200-, 
80-, and 10-micron at 
weight percent levels 
of 56/13/12/5, respec- 
tively, are shown in 
Fig. 8. We would like 
to emphasize again thet 
no manipulation of 
parameters was required; PRESSURE, PSIA 


SHENGUIVC .OP RCSL UE UEC Fig. 7. Propellant with HTPB-Type Binder 


pape Neo eran with 86% Solids and Trimodal 
Oxidizer 


current computer program. 


m— EXPERIMENTAL 
wee THEORETICAL 
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BURN RATE 
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am—— EXPERIMENTAL 
—— THEORETICAL 


BURN RATE, 


PRESSURE, PSIA 


Fig. 8. Propellant with HTPB-Type Binder 
with 86% Solids and Trimodal 
Oxidizer 


FUTURE PLANS 


In the ensuing months, we plan to further study the possibilities of 
a two surface temperature model, add a treatment for the catalyst, and 
continue studying the effects of particle size. We will also establish 
possible maximum and minimum burn rates and continue to make improvements 
in the computer program. 
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